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Abstract 

^ ■ As the Cosmic Microwave Background (CMB) radiation is observed to higher and higher 

, angular resolution the size of the resulting datasets becomes a serious constraint on their analysis. 

In particular current algorithms to determine the location of, and curvature at, the peak of 
the power spectrum likelihood function from a general iVp-pixel CMB sky map scale as 0{Np). 
' Moreover the current best algorithm — the quadratic estimator — is a Newton- Raphson iterative 

OO ! scheme and so requires a 'sufficiently good' starting point to guarantee convergence to the true 

0^ ' maximum. Here we present an algorithm to calculate bounds on the likelihood function at 

fH ! any point in parameter space using Gaussian quadrature and show that, judiciously applied, it 

Qh' scales as only 0{Np' ). Although it provides no direct curvature information we show how this 

approach is well-suited both to estimating cosmological parameters directly and to providing a 
^ ' coarse map of the power spectrum likelihood function from which to select the starting point 

, for more refined techniques. 



PACS numbers: 98.80, 98.70.V, 02.70 



1 Introduction 



Planned observations of the Cosmic Microwave Background (CMB) will have sufficient angular 
resolution to probe the CMB power spectrum up to multipoles / ~ 1000 or more (for a general 
review of forthcoming observations see Q). If we are able to extract the multipole amplitudes 
Ci from the data sufficiently accurately we will be able to obtain the values of the fundamental 
cosmological parameters to unprecedented accuracy. The CMB will then have lived up to its 
promise of being the most powerful discriminant between cosmological models |||, |3|, |^. 

Extracting the power spectrum is conceptually simple — the raw data is cleaned and converted 
into a time-ordered dataset. This is then converted to a sky temperature map, which in turn 
is analysed to find the location of, and curvature at, the maximum of the likelihood function of 
the power spectrum. In practice as the size of the dataset increases the problem rapidly becomes 
intractable by conventional methods. This is particularly true of the final step — the likelihood 
analysis of any reasonably general sky temperature map. 

An observation of the CMB contains both signal and noise 

Ai = Si + Ui (1) 

at each of Np pixels. For independent, zero-mean, signal and noise the covariance matrix of the 
data 

M = (aA^\ = (ss^\ + (nn^\ = S + iV (2) 



is symmetric, positive definite and dense. For any theoretical power spectrum Ci we can construct 
the corresponding signal covariance matrix S{Ci); knowing the noise covariance matrix N for the 
experiment we now know the observation covariance matrix for that power spectrum M{C{). The 
probability of the observation given the assumed power spectrum is then 

PiMCi) = -7^ (3) 

Assuming a uniform prior, so that 

P{Ci I A) oc P(A I Ci) (4) 

we can restrict our attention to evaluating the right hand side of equation (|3|). Unfortunately 
unless the noise covariance matrix is unrealistically simple (eg. diagonal) both direct evaluation 
1^, § and quadratic estimation ||^, |8| of the likelihood function scale as at best 0{Np) making 
them impractical for the forthcoming 10^ - 10^ pixel datasets. 



2 Bounding The Likelihood Function 

Instead of the expensive exact evaluation of the likelihood function, here we implement a much 
cheaper bounding algorithm due to Golub et al |ll|]. This method determines bounds 

C <u^f{A)u<U (5) 

for an n-vector u, symmetric positive definite n x n matrix A and smooth function / defined on 
the spectrum of A. The underlying idea is to rewrite the problem as a Riemann-Stieltjes integral 
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which is then approximated using Gaussian quadrature. Since A is symmetric it can be written in 
eigendecomposition as 

A = Q^AQ (6) 

where Q is the orthogonal matrix of normahsed eigenvectors and A is the diagonal matrix of 
increasing eigenvalues An = Aj. Writing u = Qu we have 



u^/(^)u = u^/(A)Q 

n 



„~.2 



1=1 

/(A)d/.(A)^/(/) (7) 



where the measure is a piecewise constant function defined by 

f A < Ai 

K^) = { ELi^? A, <A<A,+i (8) 

The integral /(/) can now be bounded above and below using Gauss-Radau quadrature |12] 

m 

I{f) = Y,u;if{9,) + uif{Xi)+Ri 
1=1 

m 

I{f) = Y,u:if{e{) + Vnf{K)+Rn (9) 
i=l 

with weights and vi^n, nodes Oi and Ai,n, with opposite-signed remainders -Ri,n- Given the 
resulting bounds on /(/) we can increase the number of nodes m until some convergence criterion, 
such as a maximum relative error, 

is met. Calculating such bounds scales as 0{m'n?) due to the O(n^) matrix-vector multiplication 



in using the Lanczos algorithm to calculate each of the m nodes and their weights. 
Rewriting the likelihood function of equation (^) as 

P(A I Q) oc exp (^-^ A + Tr [InM])^ (11) 

leaves A^M"^A and Tr [InM] to be evaluated. The former is already in the required form and 
can be bounded immediately. For the latter we note that 

v^/(^)v)=Tr[/(A)] (12) 

where v is a random vector whose elements are ±1 with equal probability. Generating r realisations 
of V we can calculate the bounds 

Ci<vf lnMYi<Ui l<i<r (13) 
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for each, from which we want to derive bounds on the expectation value 

^v^ InMv^ = Tr[lnM] (14) 

For each reahsation the estimator 

Xi = l {Ui + Ci)=n + 5i + ei (15) 

the sum of the expectation value )U, a random sample error and a systematic bound- width error 
ej, where 

kil < 2 (^i - A) = (16) 

(note that this is an absolute, not relative, error measure). Given the variance of the systematic-free 
data 



i=l 

= ^^ti^i-~^f (17) 

i=l 

and assuming that r is large enough for the central limit theorem to apply, our estimator X has 
student's t-distribution with (r — 1) degrees of freedom, and we can take bounds 

X - Tr-l,a S - a < < X + Tr-l,a S + 0, (18) 

with a confidence. Although the systematics prevent us from determining S itself we can calculate 
a stochastically larger quantity as follows: taking the sample variance of the midpoints 

S\X) = -l-j2{l^ + 6i + ei-fx-S-ef 
1=1 



> S^-2S 
where we have used the fact that 
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TE«f (19) 

i=l 



r 



E(^-^")'<E(e^)'<E«? (20) 



i=l i=\ i=\ 



Thus, defining 
we have 

and hence the bound 



A = 



i 



(21) 

i=i 



S'^ -2AS-S'^{X)<0 (22) 



S<A + ^A^ + 52 (X) (23) 

If the number of terms m required to evaluate each of the r bound pairs is approximately 
constant, we have a method to bound the likelihood function which scales as 0{rmn'^) overall. 
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3 Results 



The viability of this approach depends crucially on the way that the number of nodes in the Gauss- 
Radau quadrature (m) and the number of trace estimates (r) depend on the size of the dataset (n) 
and the required tightness of the bounds (e). To examine these dependencies we have applied the 
algorithm to subsets of the COBE data and a standard CDM target power spectrum. 

The number of nodes required to achieve a given accuracy clearly depends on that accuracy. 
At the extrema 

e = 1 — ^ m = 1 

e = m = n (24) 

Figure 1 shows the dependence of the number of nodes on size of the dataset for typical useful values 

e = 10^^, 10"^, 10^'^, 10^^, and 10^^. The points are from numerical experiments, evaluating 
bounds on v-^ In M v and averaging over 100 realisations of v. The solid lines are power law fits 
m (X (giving overall scaling as 0(n^+^)) with 

r 1/3 for 10-1 > e > lO'^ 
^" \ 1/2 for 10-3 > e> 10-5 ^^^^ 

Figure 2 shows the normalised 99% confidence bounds achieved on Tr [In AI] for a 1000 pixel 
dataset as the number of estimates r increases, with covergence set at e = 10"^, 10"^, 10"^, 10"^ 
and 10~^. Not surprisingly, the 10% bounds on each estimate give a rather poor overall constraint. 
However with the 1% and tighter bounds wc can determine the logarithm of the likelihood function 
to 2 - 3% with 99% confidence in as few as 20 realisations. Note that below the 1% level the bounds 
are dominated by sample error, and become essentially independent of the bound width. Figure 3 
shows the 99% confidence bounds achieved after 20 realisations at e = 10~^, showing no systematic 
variation as the size of the dataset increases. 



4 Conclusions 

We have presented an algorithm to calculate probabilistic bounds on the power spectrum likeli- 
hood function from an ATp-pixel CMB map using Gaussian quadrature which scales as between 

0(ArJ/^) and 0{N^^^) — a very significant advance on existing algorithms for the exact calculation 
which scale as 0{Np). Since lowering the convergence constraint below the 1% level gains us only 
marginally tighter final bounds at the expense of increasing the scaling power, it is not recom- 
mended for the forthcoming 10^ - 10^ pixel GMB maps. Our final algorithm of choice therefore 
gives better than 3% bounds on the logarithm of the likelihood function with 0{n'1^^) operations 
with 99% confidence. 

Since this algorithm gives no information about the local curvature of the likelihood function 
it is not as well suited as quadratic estimator techniques for searching a large multi-dimensional 
parameter space for its likelihood maximum. However for direct estimation of a small set of 
cosmological parameters this technique is certainly viable and fast. Moreover, even when the 
parameters are taken to be the multipole moments (individually or in bins), quadratic estimation, 
being a Newton-Raphson iteration, requires a starting point 'sufficiently close' to the maximum to 
guarantee convergence: the algorithm presented here is then well suited to provide a coarse overall 
mapping of the likelihood function from which to select a starting point for more refined techniques. 
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Figure Captions 



Figure 1 

A log-log plot of the scaling in the number of quadrature nodes m required to achieve bounds 
with a given relative error convergence criterion e with the size of the CMB dataset n. The points 
are obtained from numerical experiments; the lines are power law fits m cc n^. From the bottom 
of the figure reading upwards the relative errors are e = 10-^10-^10-^10-^ and 10-^ with 
corresponding power laws /3 = 1/3, 1/3, 1/2, 1/2 and 1/2. 

Figure 2 

A plot of the normalised 99% confidence upper and lower bounds achieved on Tr [In M] for a 1000 
pixel dataset against the estimator sample size r. Prom the outer limits reading inwards the line- 
pairs correspond to the relative error covergence criterion being set at e = 10~^, 10~^, 10~^, 10~^ 
and 10-^ 



Figure 3 

A plot of the variation in the normalised 99% confidence upper and lower bounds achieved on 
Tr [In M] with the relative error convergence criterion set at e = 10~^ after r = 20 estimates as the 
size of the CMB dataset n increases. 
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